library(dplyr)
library(knitr)
library(maptools)
library(rgdal)
library(TreeSegmentation)
library(sp)
library(ggplot2)
library(rgl)
library(lidR)
library(kableExtra)
knit_hooks$set(webgl = hook_webgl)
opts_chunk$set(warning=F,message=F)
#set color ramp for treeID
col = pastel.colors(200)
#set data paths
path_to_tiles="../data/2017/Lidar/"
#set cores
cores<-3
#cores<-15
shps<-list.files("../data/ITCs/test/",pattern=".shp",full.names = T)
itcs<-lapply(shps,readOGR,verbose=F)
names(itcs)<-sapply(itcs,function(x){
id<-unique(x$Plot_ID)
return(id)
})
ground_truth<-itcs[[3]]
fname<-get_tile_filname(ground_truth)
inpath<-paste("../data/2017/Lidar/",fname,sep="")
tile<-readLAS(inpath)
tile@crs<-CRS("+init=epsg:32617")
plot(tile)
You must enable Javascript to view this page properly.
plot(extent(tile),col='red')
plot(extent(ground_truth),col='blue',add=T)
silva<-silva2016(path=inpath,output="all")
dalponte<-dalponte2016(path=inpath,output="all")
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
## user system elapsed
## 1.820 0.042 1.906
## [1] "Creating tree polygons"
li<-li2012(path=inpath,output="all")
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
## user system elapsed
## 0.092 0.000 0.031
## [1] "Creating tree polygons"
watershed_result<-watershed(path=inpath,output="all")
## [1] "Computing Ground Model"
## [1] "Computing Canopy Model"
## [1] "Clustering Trees"
## user system elapsed
## 2.566 0.045 2.640
## [1] "Creating tree polygons"
plot(silva$tile,color="treeID",col=col)
You must enable Javascript to view this page properly.
plot(silva$convex)
plot(ground_truth,col='red')
plot(silva$convex,add=T)
#plot(dalponte$convex,add=T)
chm=canopy_model(silva$tile)
plot(chm,ext=extent(ground_truth))
plot(ground_truth,add=T,col='red')
plot(silva$convex,add=T)
Okay that’s not great, but let’s keep going for the moment.
Silva v Dalponte
plot(chm,ext=extent(ground_truth))
plot(silva$convex,add=T)
plot(dalponte$convex,add=T,col=rgb(0,0,255,20,maxColorValue=255))
Li versus watershed
plot(chm,ext=extent(ground_truth))
plot(li$convex,add=T)
plot(watershed_result$convex,add=T,col=rgb(0,0,255,20,maxColorValue=255))
How many tree predictions?
ptlist<-list(silva=silva$tile,dalponte=dalponte$tile,li=li$tile,watershed=watershed_result$tile)
unique_total<-sapply(ptlist,function(x) length(unique(x@data$treeID)))
df<-data.frame(Algorthm=c(names(ptlist)),Total_Trees=as.numeric(unique_total))
df %>% kable() %>% kable_styling()
| Algorthm | Total_Trees |
|---|---|
| silva | 67 |
| dalponte | 67 |
| li | 160 |
| watershed | 41 |
Each tree is assigned based on the maximum overlap. Pairwise membership is done using a Hungarian Algorithm. See clue::solve_LSAP.
assignment<-assign_trees(ground_truth,prediction=silva$convex)
#loop through assignments and get jaccard statistic for each assignment pair
statdf<-calc_jaccard(assignment=assignment,ground_truth = ground_truth,prediction=silva$convex)
ggplot(statdf) + geom_histogram(aes(IoU)) + labs(x="Intersection over union") + theme_bw()
mean(statdf$IoU)
## [1] 0.2088944
median(statdf$IoU)
## [1] 0.119387
results<-evaluate(ground_truth=itcs[[2]],algorithm = c("silva"),path_to_tiles=path_to_tiles,plot_results=T,basemap="../data/2017/Camera/L3/")
ggplot(results,aes(y=IoU,x=Method)) + geom_boxplot() + theme_bw()
results %>% group_by(Method) %>% summarize(mean=mean(IoU),median=median(IoU))
system.time(results_all<-evaluate_all(itcs=itcs,algorithm = c("silva"),path_to_tiles=path_to_tiles,cores=cores,extra=F,plot_results = T,basemap="../data/2017/Camera/L3/"))
## [1] "ITC has no overlap with cropped tile"
## user system elapsed
## 0.061 0.011 44.818
#plot results
ggplot(results_all,aes(y=IoU,x=Method)) + geom_boxplot() + theme_bw()
#Compute test statistics
test_statistic<-results_all %>% group_by(Method) %>% summarize(mean=mean(IoU),median=median(IoU)) %>% mutate(Date=format(Sys.time(), "%m/%d/%y %H:%M:%S"))
test_statistic
## # A tibble: 1 x 4
## Method mean median Date
## <chr> <dbl> <dbl> <chr>
## 1 silva 0.216 0.207 07/24/18 12:57:41
write.table(test_statistic,"Results/results.csv",append = T,col.names = F,sep=",",row.names = F)